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Abstract 

We present a new energy-stable open boundary condition, and an associated numerical algorithm, for 
simulating incompressible flows with outflow/open boundaries. This open boundary condition ensures 
the energy stability of the system, even when strong vortices or backflows occur at the outflow boundary. 

Under certain situations it can be reduced to a form that can be analogized to the usual convective 
boundary condition. One prominent feature of this boundary condition is that it provides a control 
over the velocity on the outflow/open boundary. This is not available with the other energy-stable 
open boundary conditions from previous works. Our numerical algorithm treats the proposed open 
boundary condition based on a rotational velocity-correction type strategy. It gives rise to a Robin- 
type condition for the discrete pressure and a Robin-type condition for the discrete velocity on the 
outflow/open boundary, respectively at the pressure and the velocity sub-steps. We present extensive 
numerical experiments on a canonical wake flow and a jet flow in open domain to test the effectiveness 
and performance of the method developed herein. Simulation results are compared with the experimental 
data as well as with other previous simulations to demonstrate the accuracy of the current method. Long¬ 
time simulations are performed for a range of Reynolds numbers, at which strong vortices and backflows 
occur at the outflow/open boundaries. The results show that our method is effective in overcoming the 
backflow instability, and that it allows for the vortices to discharge from the domain in a fairly natural 
fashion even at high Reynolds numbers. 

Keywords: backflow instability; energy stability; outflow; open boundary condition; outflow boundary 
condition; velocity correction 


1 Introduction 


The current work focuses on the outflow/open boundary in incompressible flow simulations and the issue of 
backflow instability, which refers to the commonly-encountered numerical instability associated 
vortices or backflows at the outflow or open boundaries. Extending our efforts on this problem 
we strive to develop effective and efficient techniques to overcome the backflow instability. 

A large class of flow problems involve physically-unbounded domains, such as jets, wakes, boundary layers, 
and other spatially-developing flows. When numerically simulating such problems, one will need to truncate 
the domain artificially to a finite size and impose some open (or outflow) boundary condition (OBC) on 
the artificial boundary. Open boundary conditions are among the most difficult and least understood issues 
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in incompressible flow simulations 


28 


63|, and have commanded a sustained interest of the community for 


decades. Among the larg e volume of works accumulated so far on this problem, the traction-free condition 
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46l | and the convective (or radiation) boundary condition [651 . l56 


are two of the more commonly used. We refer the reader to 

[^, 37. 38, 34, 23, 33, 


and to 


29 


42 


55 


22 , 


67 


62|,l8| 


631 for a review of this field up to the inid-1990s. 


for a number of other methods developed by different researchers. 


Backflow instability is a commonly encountered issue with outflows or open boundaries at moderate 
and high Reynolds numbers. Simulations have been observed to instantly blow up when strong vortices or 


backflows occur at the outflow/open boundary 


3C5 


3 3- As pointed out in 


18|, a certain amount of 


backflow at the outflow boundary appears harmless at low Reynolds numbers, but when the Reynolds number 
increases beyond some moderate value, typically several hundred to a thousand depending on the geometry 
(e.g. between Re = 300 ^ 400 for the flow past a square cylinder in two dimensions), this instability becomes 
a severe issue to numerical simulations. Commonly-used tricks for flow simulations such as increasing the 
grid resolution or decreasing the time step size are observed to not help with this instability Q Q. 

For certain flow problems (e.g. bluff-body wakes) one way to circumvent this difficulty is to employ a 
large computational domain and to place the outflow/open boundary far downstream. The idea is to allow 
for the vortices generated upstream to sufficiently dissipate before reaching the outflow boundary. This is 
feasible and computationally m anag eable at moderate Reynolds numbers. But this strategy does not scale 


with the Reynolds number 


Q,I3. 


because the domain size essential for numerical stability grows with 
increasing Reynolds number. As the Reynolds number becomes large, the needed domain size for stability 
can become very substantial. For example, in the three-dimensional direct numerical simulation of the flow 

a domain size with a wake region 50 times the 


past a circular cylinder at Reynolds number Re = 10000 [ 3 | - 
cylinder diameter in length has been used. Such a large wake region is essential for numerical stability for 
that Reynolds number, even though the far wake (beyond about 10 times the cylinder diameter) is of little 
or no physical interest and the meshes/computations in that far region are essentially wasted. 

A far more attractive approach is to devise effective open/outflow boundary conditions to overcome the 
backflow instability. Several such boundary conditions have been studied in the literature. By consideringthe 
weak form of the Navier-Stokes equation and symmetrization of the nonlinear term, Bruneau & Fabrie [y, Q] 
proposed to modify the traction condition by a term 4(n • u)^u, where u and n are respectively the velocity 


and the outward-pointing unit vector normal to the outflow boundary, and (n • u) is defined as n • u if 
n • u < 0 and is set to zero otherwise. We refer to e.g. [^, 3 applications of this boundary condition 


in later works. A traction condition containin g a term (n • u) u, which is very similar to that of 


without the 4 factor, has been employed in [l 

has also been considered in 0 . By considering the energy balance of the system, we have proposed 


0,3 


but 


50 . 59 . 27 . 361. Note that a form /3(n • u) u where 0 < /3 < 1 


in 0| 


2 
















































































a boundary condition involving a term i|upn0o(n, u), where |u| is the velocity magnitude and 0o(n, u) is 
a smoothed step function about n • u (see Section [2] for definition). While the role of the term 0o(n, u) can 

e OBC in [l^ is very different from 


those involving (n • u)u of the previous works 


in [ 4 1 employs a penalization of the tangential velocity derivative to allow for improved energy balance. Very 


E, 1,1, M, 27, 3- 


Another boundary condition developed 


recently we have proposed in 


18| a family of open boundary conditions, having the characteristic that they 


all ensure the energy stability of the system even in situations where strong vortices of backflows occur at the 


E, 1, 27, 36, 14| 


as particular 


outflow/open boundary. This family of boundary conditions contains those of 
cases, and more importantly provides other new forms of open boundary conditions. Several of those forms 
have been investigated in relative detail in [l^. 

It is observed that, while some of the above open boundary conditions have existed in the literature for 
some time, their adoption in production flow simulations appears still quite limited. This is perhaps in part 
owing to the challenge associated with the numerical implementation of these boundary conditions. All the 
aforementioned boundary conditions for tackling the backflow instability couple together the velocity and 
the pressure, and it is not immediately clear how to implement them in numerical simulations. This seems 
to be exacerbated by the fact that, when these boundary conditions are originally proposed, for most of 
them their numerical treatments are not discussed or not adequately discussed, especially in the context of 
the commonly-used splitting or fractional-step type schemes for incompressible flow simulations. It is noted 


that in the more recent works 
type strategy 
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two splitting-type schemes, respectively based on a velocity-correction 




are presented to deal with the energy-stable 


and a pressure-correction type strategy 
open boundary conditions developed therein. These algorithms de-couple the computations for the pressure 
and the velocity in the presence of open/outflow boundaries. 

The objective of the current paper is twofold. First, we present a new energy-stable open boundary 
condition that is effective in overcoming the backflow instability for incompressible flow simulations. This 
boundary condition involves an inertia (velocity time-derivative) term, and can be shown to ensure the energy 
stability of the system even in the presence of backflows or vortices at the open/outflow boundary. It does 


Q- 


not belong to the family of open boundary conditions discussed in lia . If no backflow occurs at the outflow 
boundary, this boundary condition can be reduced to a form that can be analogized to the usual convective 
boundary condition. Hence, we refer to it as the convective-like energy-stable open boundary condition. 
The current open boundary condition has a prominent feature: it provides a control over the velocity at the 


open/outflow boundary. In contrast, the family of energy-stable open boundary conditions from 


181 and 


the other aforementioned boundary conditions to address the backflow instability do not provide any control 
over the velocity at the open/outflow boundary. Therefore, as the vortices pass through the outflow/open 
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boundary, the current boundary condition can lead to smoother velocity patterns in regions at or near the 


outflow boundary when compared to that of 


m- 


Second, we present an efficient numerical algorithm for treating the proposed open boundary condition. 
Our algorithm overall is based on a rotational velocity-correction type splitting approach, and the key issue 
lies in the numerical treatment of the inertia term in the open boundary condition. At the pressure sub-step 
our scheme leads to a Robin-type condition for the discrete pressure on the outflow boundary, and at the 
velocity sub-step it leads to a Robin-type condition for the discrete velocity on the outflow boundary. In 


Ei, 18| 


contrast, the algorithms of |14l. Il8j both impose a pressure Dirichlet type condition on the outflow boundary 
at the pressure sub-step and a velocity Neumann type condition on the outflow boundary at the velocity 
sub-step. The current algorithm is simpler to implement with spectral-element (and also finite-element) type 
spatial discretizations, because there is no need for the projection of pressure Dirichlet data on the outflow 
boundary as required by the algorithms of [l4 , 1^ . 


We would like to point out that, by using an idea analogous to that of m , one can generalize the current 
open boundary condition to a family of convective-like energy-stable open boundary conditions, which will 
provide other new forms of OBCs; see the discussions in Section ETT] in this regard. The numerical algorithm 
presented herein with no change can be applied together with this family of convective-like energy-stable 
OBCs. 

The novelties of this work lie in two aspects: (i) the convective-like energy-stable open boundary condition, 
and (ii) the numerical algorithm for treating the proposed open boundary condition. The rotational velocity- 
correction scheme for discretizing the Navier-Stokes equations employed here has also subtle differences than 


Q 


that of M in the numerical approximations of various terms, although both can be classified as velocity- 
correction type schemes. 

The open boundary condition and the numerical algorithm developed herein have been implemented and 


tested using high-order C° spectral elements for spatial discretizations (6J, l41 , 


75j | . The implementations 


discussed in the paper without change can also be used for finite element discretizations. It should be noted 
that the open boundary condition and the numerical algorithm are very general and are not confined to a 
particular spatial discretization technique. They can also be implemented using other spatial discretization 
methods. 
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2 Convective-Like Energy-Stable OBC and Algorithm 

2.1 Convective-Like Energy-Stable Open Boundary Condition 

Consider a domain 12 in two or three dimensions, and an incompressible flow contained within this do¬ 
main. Let 912 denote the boundary of the domain 12. This flow problem is then described by the following 
incompressible Navier-Stokes equations in non-dimensional form: 


u • Vu -1- Vp — = f, 

at 

(la) 

V • u = 0, 

(lb) 


where u(x,f) and p(x, f) are respectively the normalized velocity and pressure fields, f(x, t) is an external 
body force, and x and t are the spatial coordinate and time. The constant v denotes the normalized viscosity, 
z/ = -^, where Re is the Reynolds number defined by appropriately choosing a length scale and a velocity 
scale. 

We assume that the domain boundary 912 consists of two types: 

912 = 912d U 912o, 912d n 912o = 0. (2) 

In the above 912^ denotes the Dirichlet boundary, on which the velocity is given 


u = ’w(x, t), on 912d 


(3) 


where w is the boundary velocity. On the other hand, on 912o neither the velocity u nor the pressure p is 
known. We refer to 912o as the open (or outflow) boundary. How to deal with 912o in numerical simulations 
is the focus of the current work. 

We propose the following boundary condition for the open boundary: 


1 

i^Dq— - pn + lyn-Wu- - |u|^ n -f- (n • u) u 0o(n, u) = f{,(x, t). 


on 912o. 


(4) 


In the above equation, n is the outward-pointing unit vector normal to the boundary 912o. Dq is a chosen 
non-negative constant {Dq ^ 0), which has been normalized by ^ {Uq denoting the characteristic velocity 
scale) and is non-dimensional, ft, is a prescribed vector function on 912o for the purpose of numerical testing 
onl y, a nd it is set to ft, = 0 in actual simulations. 0o(n,u) is a smoothed step function about (n • u) given 

by y. 


^ N 1 A , n • u 
e„(n,u) = -h-la„h — 


(5) 


where i5 > 0 is a non-dimensional positive constant that is sufficiently small. As discussed in [1^ , d controls 
the sharpness of the smoothed step function, and it is sharper if <5 is smaller, and that the simulation result 
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is not sensitive to S when it is sufficiently small. As (5 —>■ 0, 0o(n,u) approaches the step function, that is, 


lim 0o(n, u) = 0so(n, u) = 


1 , if n • u < 0, 


( 6 ) 


s^o , I otherwise. 

A prominent characteristic of the open boundary condition (j4]) is the inertia term ^ involved therein. 
One can note that for Dq = 0 the inertia term vanishes and the boundary condition Q will be reduced to 
the so-called boundary condition “OBC-C” that has been studied in . In the current work we concentrate 
on the cases with Dq > 0. 

The open boundary condition (|3]), with ff, = 0 and when S is sufficiently small, ensures the energy stability 
of the system. To illustrate this point, we look into the energy balance equation for the system consisting of 
(fTal) and (fTbl) : 


^ ^ f ||Vu||^ 

^ JQ 


dt 


f • u 


/n 


/n 


' dVld 


n • T • u — - |u|^ n • u 


-b 


Idila 


n • T • u — - |u|^ n • u 


(7) 


where T = —pi -b i/Wu and I denotes the identity tensor. We assume that fh = 0 in Q and (5 —>■ 0 in 
0o(n,u). Then by employing the condition (|1]) on dHo, the last surface integral on the right hand side of 
© becomes 


idna 


n -T u - i|u|^n u 


dt. 

JdQo 

d 

r 


Jacia 


uDo^ |u|^ -b 
lyDo^ |u|^ - 


’dQo 


'dQo 


^ |u|^ (n • u) [20so(n,u) - 1] 


( 8 ) 


u n-u 


as S —^ 0. 


It then follows that the energy balance equation can be transformed into 

1 

I dQo 

f U 


dt 


^\uf + 1 ^ a 


0 


2 


= -i. / ||Vu|| 

Jn 


n-Tu—-|u|^n-u| — 


1 


(9) 


u n-u 


as 5 —>■ 0. 


Jn JdQd \ ^ / Jdfrio ^ 

Therefore, the open boundary condition given by o, when ff, = 0 and 6 is sufficiently small, ensures the 
energy stability of the system (in the absence of external forces), even if strong vortices or backflows occur 
(i.e. n • u < 0) on the open boundary dflo- Note that, because the velocity u is given on d^ld, the surface 
integral term on dlld in equation (jO)) will not pose a numerical instability issue. 

It is instructive to compare the energy balance equations for the current open boundary condition and 
for the open boundary conditions introduced in 1^. Let us assume for now that there is no external 


body force f and that u = 0 on the Dirichlet boundary dild- Then equation © implies that the sum 
(b I |u|^ -b i^Dq fg^ I will not increase over time. For Dg > 0, the energy balance relation provides 
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an upper bound for the total energy ^ |u|^ and for the quantity fg^ ^ |u|^ with the current open boundary 
condition. This provides a control over the velocity u on the outflow boundary dfio- On the other hand, 
with the open boundary conditions from 18| , the energy balance equation involves only the total energy, and 


there is no control over the velocity on the outflow boundary. This is a key difference between the current 
open boundary condition and those from Q . Thanks to this characteristic, the current open boundary 
condition can lead to qualitatively smoother velocity patterns at/near the outflow boundary as vortices pass 
through. This point will be illustrated in Section using numerical simulations. 

In addition to the open boundary condition (|4]), we will also consider the following boundary condition. 


clu 

vDq— -pn + j/n-Vu = 0, on 

ot 


or equivalently for T>o > 0 
i9u 1 (?u 1 


pn, on d^lo- 


( 10 ) 


( 11 ) 


dt Dq dn vDo 

The difference between this boundary condition and (|3]) lies in that this boundary condition does not ensure 
the energy stability when backflow occurs on the open boundary dilo- In contrast, the condition (jT]) ensures 
the energy stability even in the presence of backflows at the open boundary. Note that for Dq = 0 the 
condition will be reduced to the traction-free boundary condition. One also notes that Equation m 
is reminiscent of the usual convective boundary condition (together with p = 0), 


9u _ 

— -k = 0, on d^o 


( 12 ) 


p = 0, on dflo 

where Uc denotes a convection velocity. Because of this resemblance to the convective boundary condition 
we will refer to the boundary condition (|3]) as a convective-like energy-stable open boundary condition. 

The analogy between the current boundary condition and the usual convective boundary condition sug¬ 
gests that in the boundary condition (|T]) the parameter plays the role of a convection velocity scale at the 

velocity in the usual convective boundary condition 


outflow boundary. Different choices for the convection 
have been considered in a number of studies (see e.g. 


63|, l8|), which can provide a guide for the choice of 


Dq in the boundary condition (jT|). For a given flow problem, one can first perform a preliminary simula¬ 
tion using the boundary condition Q with Dq = 0 to obtain an estimate of the convection velocity scale 
C/c > 0 at the outflow boundary, and then carry out the actual simulation by setting Dq = in (|T]). Our 
numerical experiments in Section [3] indicate that the difference in the Dq values has very little or no effect 
on the computed global flow quantities such as the force coefficients. The main effect of Dq appears to be on 
the qualitative flow characteristics local to the outflow boundary. An improved estimate of the convection 
velocity (and hence an improved Dq value in (|3])) will allow the vortices or other flow features to discharge 
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from the domain more smoothly and in a more natural fashion. We will look into the effects of value in 
the open boundary condition Q in more detail in Section [3] 


Remarks By employing an idea similar to that of 18| , we can come up with other forms of convective- 
like energy-stable open boundary conditions. We briefly mention several of them below: 


d\i 

vDq— — pn + vn- Vu — [(n • u) u] 0o(n, u) = 0, on 9f2o; 


(13a) 


7^ (?U 

vDq— - pn -I- • Vu — 

ot 


lul^n 


0 o(n, u) = 0, on 


(13b) 


9u 

vDq— -pn -I- • Vu — 

ot 


- (n • u) u 


0o(n,u)=O, on cirio; 


(13c) 


9u 

vDq— - pn + vn- Vu — 

at 


luP n 


0o(n,u)=O, on Sflo; 


(13d) 


^U. 1 

uDq— - pn + z/n-Vu - |uPn+(n-u)u 0o(n,u) = O, on d^lo- 

ot 4 L 


(13e) 


We would also like to mention the following more general form (analogous to [l^), which contains the 
boundary condition (j4|) and those represented by (jl3a|l - (jl3ell as particular cases, 


(?U 

vUq— -pn -I- • Vu — 

at 


[0 + a 2 ) ^ |u|^ n -f (1 - 6» -H ui) i (n • u) u 


where 0, ai and a 2 are constants satisfying the conditions 

0 ^ 0 ^ 1, ai > 0, 02 ^ 0. 


0o(n,u)=O, on Srio, (14) 


(15) 


Note that the general form m ensures the energy stability of the system as (5 —?► 0. In this case the energy 
balance relation is given by the following expression, 




1 


afla 


= -^ / l|Vu| 
Jn 


P+ / f u 

Jn 


' dCld 


n • T • u — i |u|^ n • u 


+ [ i |u|^ (n • u) [(1 + ai+ 02 ) 0so(n ■ u) - 1] 


(16) 


laQo 

/ ||Vu| 
Jn 


P+ / f u- 

Jn 


land 


n • T • u — i |u|^ n • u ) , as i5 —>■ 0. 


Apart from the boundary conditions discussed above, we impose the following initial condition for the 
velocity, 


u(x,t = 0) = Ui„(x), 


(17) 















where Ui„ is the initial velocity field satisfying equation (llbl) and compatible with the boundary condition 
(l3|) on dfld at t = 0. 

2.2 Algorithm Formulation 

The equations (El) and (llbl) . the boundary condition ([3]) on dfld, and the boundary condition (jd]) on dflo, 
as well as the initial condition ED, together constitute the system that need to be solved in numerical 
simulations. 

We next present an algorithm for numerically simulating this system, with emphasis on the numerical 
treatment of the open boundary condition (j3]) . Let 


1 r 2 

E(n,u) = -yu| n+(n-u)u 0o(n, u), 


and we re-write equation (1T|) into a more compact form, 

dll 

i'Dq— -pn + z/n • Vu — E(n, u) = f^. 

ot 


(18) 


(19) 


We will concentrate on the algorithm and implementation for Dq > 0 in m in this and the next subsections. 
In section 12.41 we will briefly discuss how to deal with the case Dq = 0, when the current open boundary 
condition is reduced to a form already studied in [l^ . 

Let n ^ 0 denote the time step index, and (•)" denote (•) at time step n. Define u° = Ui„. Then, given 
u" we compute ^ u"'+^) in a de-coupled fashion as follows: 

For p"+h 


{^n+1 _ 

2E__ -+ + Vp"+^ -b X V X u*’"+i = r+1 


(20a) 


V • u"+i = 0 


(20b) 


n • = n • on d^d 


(20c) 


Th~\~ 1 *** 

vDq——— -- n — -b vn ■ • n — n • E(n, u*’"^^) = • n, on 9Do 


At 


For u"+^: 


At 


(20d) 


(21a) 


u"+i=w"+\ onSDd 


(21b) 
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vDi 


7 ou"+^ — u 

' aT 


-P 


n+1 


n + i^n-Vu”+i-E(n,u*'"+^) + i/(V-u*’”+i)n = 4 "+\ on d^o- ( 21 c) 


In the above equations, At is the time step size, n is the outward-pointing unit vector normal to the boundary, 
and is an auxiliary variable approximating Let J (J = 1 or 2) denote the temporal order of 

accuracy of the algorithm. Then is a J-th order explicit approximation of given by 


*,n+l _ 


u", 

2u” - u"-h 


J=l, 

J = 2. 


( 22 ) 


The expressions i( 7 ou”+^—u) and i( 7 ou"+^ —u) are approximations of by a J-th order backward 


differentiation formula, and u and 70 are given by 

7o = 


2u” - 


J=l, 
J = 2, 


J = l, 
J = 2. 


(23) 


Note that E(n, u) is given by (ITS)) . 

One can observe that the overall structure of the above algorithm represents a rotational velocity- 
correction type strategy (see [^, 16, 17, 141) for de-coupling the computations of the pressure and velocity. 
While both belong to velocity correction-type schemes, the scheme here is somewhat different from the one of 
[ 1 ^ . Note that in the pressure sub-step we have approximated all terms at time step (n-l-1) with the current 
scheme Isee eauation (I20al)l. In contrast, in D the viscous and the nonlinear terms are approximated at 
time step n rather than at (n -I- 1) in the pressure sub-step, and correspondingly some correction terms are 
incorporated in the subsequent velocity sub-step. The current treatment of various terms is observed to yield 


smaller pressure errors and comparable velocity errors compared to that of 


Q- 


The inertia term vDq ^ in the boundary condition (I19|) demands some care in the temporal discretization. 
The discrete equation (I20dl) in the pressure sub-step stems from an inner product between the directional 
vector n and the open boundary condition (IT^ on the outflow boundary dflg- Note that the ^ term and 
the pressure term have been treated implicitly in (I20d|) . and in particular that ^ is approximated using 
u”+^ (instead of u”+^) here in this discrete equation. This point is crucial, and it effectively leads to a 
Robin-type condition for the pressure on d^lo because of the equation (I20al) . An explicit treatment of 
the ^ term in (I20dl) would seem more attractive and would result in a Dirichlet type condition for 
on drio, just like in the scheme of [l^. This treatment however does not work, and is unstable unless Dq 
is very small. At the velocity sub-step, in the discrete equation (l21cD on d^o we have treated the terms ^ 
and n • Vu implicitly, and note that ^ is approximated using here. These numerical treatments give 
rise to a Robin-type condition for the discrete velocity on dflo, noting that in (I21cl) is already 

explicitly known from the pressure sub-step. Note also that in the discrete equation ()21c|l an extra term 
j^(V • u)n has been incorporated in the formulation. 
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We would like to point out that the algorithmic formulation given by (I20al) - (l21cl) can be used together 
with the general form of convective-like energy-stable open boundary condition (|14l) . by setting E(n, u) in 
the algorithm as follows 


E(n,u) = 


(6» -f 02 ) ^ |u|^ n -f (1 - 6» -f oi) i (n • u) u 


0o(n,u). 


(24) 


2.3 Implementation with CO 

Spectral Elements 


We employ C^-continuous high-order spectral elements (6J, 1^ 


75l| for spatial discretizations in the current 


work. Let us next look into how to implement the algorithm, given by (I20al) - (l21cl) . using C° spectral elements. 
The discussions in this subsection with no change also apply to (7° finite element implementations. 


As noted in several previous works 


Q, 17 


id. Il4l lll|. the complication in the implementation with 
elements stems from the high-order derivative terms such as V x V x u involved in this type of algorithm, 
because such terms cannot be directly computed in the discrete function space of (7° elements. We can 
eliminate such complications by looking into the weak forms of the algorithm. In addition, we will eliminate 
the auxiliary velocity from the final form of the algorithm. 

We first formulate the weak forms of the algorithm in the spatially continuous space. Let g(x) denote a 
test function. By taking the inner product between Vg and equation (I20a|) and integrating by part, we 
have 


• Vg 


At. 


n • u 


n+1 _ 


q= / G 


n+1 


Vq-v 


dQ.o 


In 


IdndUdQo 


n X • Vg 


At JgQj 


(25) 


n • w"+^g. 


Vg, 


where u; = V x u is the vorticity. 


QK+l _ £n+l _j_ ^ ^ '^^*•'"'+1 

At 


(26) 


and we have used equations (I20bl) and (I20c|) . the divergence theorem, and the identify 


/ V X o) • Vg = / V • (cp X Vg) = / n x cp • Vg. 

In Jn Jan 


(27) 


According to equation (I20dl) . n • u"+^ can be expressed in terms of and other explicit quantities on 
dilo- We therefore can transform equation (l25l) into the final weak form for the pressure 

1 


/ V++' • Vg , 

In Jan. 


p^+\ = / G"+^ -Vq-iy n X ■ Vg 

Jn Janduana 


Jana 

At, 


-^n • u -f [z/n • VW’^+i • n - n • E(n, u*’”+i) - 1+^ ' n] | g (28) 


At 

n • w^+^g, Vg. 


and 


11 


























We next sum up equations (I21a|) and (I20al) to obtain 


- V2u"+^ = - (G’"+i - Vp”+M . (29) 

vAt V ' ' 

Let </3(x) denote a test function that vanishes on d^d- Taking the inner product between (p and equation 
(|29ll and integrating by part lead to 


70 / 

Jq 


u"+V - 


V(p ■ Vu^+i - 


[ n • Vu"+V = - / (G"+i - Vp"+i) ip, 

Jdno 


\/p, 


(30) 


where we have used the divergence theorem, and the fact that n-= 0 thanks to the requirement 
that = 0 on dfid- According to equation (I21cll . n • can be expressed in terms of u"+^ and other 

explicit quantities on dfio- We therefore can transform (I30|) into the final weak form for 

Vp ■ Vu”+i 


70 

vAt 


[ u"+V+ / Vv9-Vu”+1 + ^^ / u"+V=- f (G"+1 - Vp”+^) 

Jn Jq. JdQo ^Jq 

+ [ (?!"+- + E(n> ^*’'‘+1) + - ly {V ■ n] ] p, Vp. 

JdQo [Atu ) 


(31) 


The weak forms of the algorithm in the continuum space consist of equations (1^51) and m, together 
with the velocity Dirichlet condition (I21b|) on d^ld- The auxiliary variable does not appear in the weak 
form and is not explicitly computed. These equations in weak forms can be discretized using C° spectral 
elements (or finite elements) in a straightforward fashion. 

Let flfi denote the domain discretized using a spectral element mesh, and = dfldh U dfloh denote 
the boundary of flh, where dfldh and dfloh respectively represent the discretized dfld and d^lo- Let X?i C 
(where d = 2 or 3 is the spatial dimension) denote the approximation space of the discretized 
velocity and define Xho = { u G H^{Qh) ■ v\gaah = 3 } . Let Mh C denote the approximation 

space of the discretized pressure We take the test function q of equation (1^ from Mh, and take the 

test function p of equation (I^Tl) from Xho- In the following let {■)h denote the discretized version of the 
variable (•). Then the discretized version of equation (1^51) is: Find G Mh such that 


[ . Vg, + ^ / pl+\h = f G^i • Xqh - ly [ 

JQh JdVtoh Jd 


*,n+l 

nh X 


Vg;^ 


+ 




JdQoh 
70 [ 

AtJoQ, 


dQdhUdQoh 

*,n+l', pn+1 


lyrih ■ Vu*’”+ -Uh-nh- E(n^, " • Uh 


qh 


T^h ■ \/qh G Mh. 


The discretized version of equations (IHTl) and (I21bl) is: Find G Xh such that 


70 

uAt 


'TO-^O f 72-1-1 1 f /^n+l 

+ —^ V’h = - [G^ - Xpt^ ) ph 

JdQoh ^ JQh 

jga ^ \ + E(n,„u;;’"+^) + 4"+^ -v{v- nh] | Ph, 


u^^^Ph + Vph- Vu 
J Qh 
Do 


' dQoh 

'iph G Xho 


(32) 


(33) 
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together with 


ur‘=wr'. 


on d^dh- 


(34) 


Our final algorithm therefore consists of the following operations within a time step: (i) Solve equation 
for (ii) Solve equation (l33)) . together with the Dirichlet condition (l34l) on dfldh, for The 

computations for the pressure and the velocity are de-coupled, and the computations for the three components 
of the velocity are also de-coupled. All the terms on the right hand sides of equations (15^ and (I55|) can 
be computed directly using (7° spectral elements. Note that the auxiliary velocity is not explicitly 

computed. 

We employ equal orders of expansion polynomials to approximate the pre ssure and the velocity in the 


n 


17 . 10 , 12 . 14, 181. Note that in 


current spectral-element implementation, similar to our previous works ml 
all the numerical simulations and flow tests of Section |3] we have used the same polynomial orders for the 
pressure and the velocity. We refer the reader to the equal-order approximations for the pressure/velocity 
by other researchers in the literature 


fl 68, 31, 41, 47, 3, 1, 3- 


2.4 The Case of Dq = 0 in Open Boundary Condition 

So far we have focused on the case Dq > 0 in the open boundary condition (1191) . In this subsection we briefly 
discuss the case Uq = 0 in the boundary condition. 

As noted in Section [24] with Dq = 0 the boundary condition (jd] is reduced to a form (so-called “OBC- 


C”) that is already studied in [3- One can therefore employ the algorithms from 


case. Note that the algorithm presented 


in [3 is wi 


Q or Q 


to treat this 


n 


with respect to the open boundary condition having a 


of open boundary conditions given in 


form corresponding to the so-called “OBC-E” in Il8| . But the algorithm of [3 also applies to other forms 


1^. 


With Dq = 0 the essential difference when compared with the scheme presented in Section 12.21 lies in 
that, in the pressure sub-step the pressure condition on the open boundary will now become of Dirichlet 
type rather than Robin type, and in the velocity sub-step the velocity condition on the open boundary will 
become of Neumann type rather than Robin type. 


We now briefly mention an algorithm for Dq = 0, as an alternative to the one presented in 
discretize the governing equations and the boundary conditions as follows: 

For 


Cl- 


We 


Use equations (I20all . (I20b|l . and 


Cl); 


n+1 _ 


= i^n ■ Vu*’"+^ • n - n • E(n, - f, 


n+1 

b 


n, 


on 


(35a) 
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For u"+^: 

Use equations (I21al) and (I21b|) : 


n • Vu”+^ = - [p”+in + E(n, u*’”+i) - {V ■ u*’"+i) n + f”+i] , on 


(36a) 


D 


The difference between this algorithm and that of M lies in that, in the pressure sub-step of this algorithm 
we have approximated all terms at the time step (n -I- 1) and in the velocity sub-step no correction terms 


are involved. On the other hand, in M certain terms are approximated at time step (n -I- 1) and the other 
terms are approximated at step n in the pressure sub-step, and in the velocity sub-step several correction 
terms are involved as a result. 

The weak forms of this algorithm can be obtained using a procedure similar to that of Q . Let iLpio(O) = 
{ w G iL^(O) : = 0 } , and q G Hpo{Q) denote the test function. Then the weak form for is 


Vp' 


n+l 


Vg = / G"+l -S/q-u 


n X u; 


>k,n+l 




Vg 


At 


n • w 


.n+l 


dQd 


Vg G 


(37) 


where is given by (1^ . Let i7^g(0) = { ?; G H^{Q) : ^ T’ ^ ^uoi^) denote the test 

function. Then the weak form for is 


70 

uAt 


u’^+V. 


Vp • Vu”+i = 


- [ (G"+i - Vp"+i) p 

^ Jn 


+ - [ [p”++ + E(n,u*’"+i)+f”+i-z/(V-u*’"+i)n] p, 

^ JdQo 


(38) 


The algorithm involves the following operations within a time step: (i) Solve equation (1571) . together 
with the pressure Dirichlet condition (I35al) on dfto, for (ii) Solve equation (1551) . together with the 

velocity Dirichlet condition (121bl) on dild, for When imposing the pressure Dirichlet condition (l35aD 

on dfto using spectral elements (or finite elements), a projection of the pressure Dirichlet data to the 
i7^(i90o) space is required because of the velocity gradient term involved in the equation; see Q for more 
detailed discussions in this regard. We have implemented the above algorithm, and the numerical experiments 
reported in Section [5] corresponding to Dq = 0 are performed using this algorithm. 


3 Representative Numerical Tests 

In this section we consider several flow problems with open/outflow boundaries and employ two-dimensional 
simulations to demonstrate the effectiveness and performance of the open boundary condition and the nu¬ 
merical algorithm developed in the previous section. At large Reynolds numbers the presence of vortices and 


14 














>> 

c 

D 

O 

2 

o 

‘C 

Q 


B 


Dirichlet boundary 


Dirichlet boundary 


Open boundary 


C3 

c 

D 

O 

c 

(U 

a. 

O 


Dirichlet boundary 


D 


(a) 



Figure 1: Spatial/temporal convergence rates: Flow configuration and boundary conditions (a), and 
and errors as a function of the element order with fixed At = 0.001 (b) and as a function of At with a 
fixed element order 16 (c). 


backflows at the open/outflow boundary makes these problems very challenging to simulate. We will look 
into the spatial and temporal convergence rates of the algorithm, and compare current simulation results 
with the experimental data as well as other simulations from the literature. The results show the effectiveness 
of the proposed method for dealing with the backflow instability. 

3.1 Convergence Rates 

In this subsection we study the spatial and temporal convergence rates of the algorithm presented in Section 
12.21 by considering an analytic solution to the Navier-Stokes equation together with the open boundary 
condition proposed in Section [2Tl 

Figure HJa) shows the problem setting. Consider the rectangular domain ABDC defined by 0 ^ a; ^ 2 
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and — 1 ^ y ^ 1, and the following analytic solution to the Navier-Stokes equations, (flal) and (fTb| . 
u = 2 cos Try sin TTXsint 

V = —2 sin Try cos nx sin t (39) 

p = 2 sin Try sin ttx cos t 

where u = {u,v). We use a characteristic velocity scale Uq = 1 and a non-dimensional viscosity it = 0.01 
for this problem. The external body force f(x, t) in (jlal) is chosen such that the expressions given by (1391) 
satisfy the equation (fTa)l . It is noted that the analytical solution (IM)) employed here has been used for the 
convergence tests in previous works [l4 , 1^ . 

To simulate the problem we discretize the domain using two equal-sized quadrilateral elements {ABFE 
and EFDC) along the x direction. On the sides BD, AB and AE we impose the Dirichlet condition (jS)), 
where the boundary velocity w(x, t) is chosen according to the analytic expressions given in ([39]). On the 
sides EC and CD we impose the open boundary condition (|4|), in which we set Dq = 1.0 and 5 = -^ and 
ft,(x, t) is chosen such that the velocity and pressure expressions given by (l3^ satisfy the condition (|4|) at 
these boundaries. 

We integrate the Navier-Stokes equations (ITa)) - (llb|) using the scheme presented in l2.2l in time from t = 0 
to t = tf (tf to he specified below). Then we compute the errors of the numerical solution at t = tf 
against the analytic expression given in (1391) . The element order or the time step size At has been varied 
systematically, and the errors are collected and monitored to study the convergence behavior of the method. 

Let us first look into the spatial convergence behavior. In this group of tests we fix the time step size 
at At = 0.001 and the integration time at tf = 0.1 (100 time steps), and then vary the element order 
systematically between 2 and 20. The numerical errors corresponding to each element order have been 
computed and monitored. Figure HKb) shows the L°° and errors of the velocity and the pressure as a 
function of the element order from these tests. As the element order increases but within order 12, all the 
numerical errors are observed to decrease exponentially. When the element order increases to 12 and beyond, 
the error curves are observed to level off at a level ^ 10“^ for this problem. The saturation of the total 
numerical error is because the temporal truncation error becomes dominant when the element order becomes 
large. These results demonstrate the spatial exponential convergence rate of our method. 

The temporal convergence behavior of the method is demonstrated by Figure [TJc), in which we plot 
the L°° and errors of the flow variables as a function of the time step size At. In this group of tests 
the integration time is fixed at tf = 0.5, the element order is fixed at 16, and At is varied systematically 
between At = 0.1 and At = 0.000390625. The convergence appears somewhat not very regular when At 
is above 0.025, especially in terms of the L°° error norms. As At decreases below 0.025, one can observe a 
second-order convergence rate in time for all the flow variables. 
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Figure 2: Contours of vorticity at Reynolds numbers Re = 30 (a) and Re = 200 (b). Dashed curves indicate 
negative vorticity values. 


The results of this section suggest that for problems involving open boundaries the method presented in 
Section [2] exhibits an exponential convergence rate in space and a second-order convergence rate in time. 

3.2 Flow Past a Circular Cylinder 

In this subsection we consider a canonical wake problem, the flow past a circular cylinder, in two dimensions 
to test the performance of our method. The goal is to demonstrate the accuracy of the method by comparison 
with the experimental data, and to demonstrate its effectiveness in dealing with the backflow instability as 
the Reynolds number becomes large. 

This flow problem has been used in [l^ to test a different set of open boundary conditions and an 
associated pressure correction-based numerical algorithm. The flow configurations employed in the current 


work largely follow those of 18|. It should be noted that the open boundary condition and the algorithm 
being tested here are very different from those of [l^ . 

Specifically, we consider a circular cylinder of diameter d, and a rectangular domain containing the 
cylinder, —5d ^ x ^ L and — lOd ^ j/ ^ lOd, where x = L is the right domain boundary to be specified 
below. The center of the cylinder coincides with the origin of the coordinate system. Four flow domains 
have been considered with different wake-region sizes. They respectively correspond to L/d = 5, 10, 15 and 


20, and are chosen in accordance with 


18l| . The flow domain with L/d = 10 is illustrated in Figure dja). 


On the top and bottom domain boundaries (y = ±10d) we assume that the flow is periodic. So the 
configuration in actuality corresponds to the flow past an infinite array of cylinders aligned in the y direction. 
On the left boundary (x = —5d) a uniform flow comes into the domain with a velocity u = (u, v) = (f7o, 0), 
where f7o = 1 is the characteristic velocity scale. The right domain boundary (x = L) is assumed to be open. 
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where the fluid can freely move out of the domain and backflow may occur depending on the flow regime 
and the domain size. 

In order to simulate the problem, we discretize the domain using a mesh of quadrilateral spectral elements. 
The meshes for these four domains respectively contain 968, 1228, 1488 and 1748 quadrilateral elements. 

In simulations we impose the periodic condition at y/d = ±10, and the velocity Dirichlet condition ([3|) 
at the inflow boundary x = —5d with a boundary velocity w = (Uo,0). On the cylinder surface a velocity 
no-slip condition is imposed, i.e. the Dirichlet condition ([J) with w = 0. At the open (outflow) boundary 
X = L we impose the open boundary condition ((4|) with fh = 0 and 6 = 

We employ the algorithm developed in Section[2]to solve the incompressible Navier-Stokes equations. All 
the length variables are normalized by the cylinder diameter d and the velocity is normalized by Uq- So the 
Reynolds number for this problem is defined by 





Upd 


(40) 


where Vf is the kinematic viscosity of the fluid. A range of Reynolds numbers (up to Re = 10000) has been 
considered. We use an element order 6 for Reynolds numbers below 100, and an element order 8 for higher 
Reynolds numbers. For selected Reynolds numbers we have also performed simulations with even larger 
element orders (up to 12), and we observe that the difference in the results when compared with the element 
order 8 is small. The non-dimensional time step size is UoAt/d = 10“^ for Reynolds numbers below 100 and 
UoAt/d = 2.5 X 10“^ for higher Reynolds numbers in the simulations. 

As discussed in Section [21 the analogy between the open boundary condition m and the convective 
boundary condition m suggests that represents a convection velocity. For the majority of simulations 
in this section we employ the average velocity at the outflow boundary, C/q, as this convection velocity and 
set Do = in the open boundary condition (|4|). This is the default Do value for the results reported in this 
section. For several selected Reynolds numbers we have also investigated the effects of Do on the simulation 
results. Results corresponding to the other Do values will be explicitly specified. 

The cylinder wake can be classified into several regimes, exhibiting a variety of flow features. These have 
been expounded in the review paper [^ . For Reynolds numbers below about Re = 47 the cylinder flow is 
two-dimensional and at a steady state. As the Reynolds number increases beyond this value, the cylinder 
wake becomes unsteady and is characterized by vortex sheddings. It remains two-dimensional for Reynolds 
numbers up to about Re = 180. As the Reynolds number increases beyond Re 180, the cylinder wake 
develops an instability and the physical flow becomes three-dimensional. More complicated flow features and 
turbulence develop in the cylinder wake when the Reynolds number increases further. In Figures |2Ka) and 
(b) we plot contours of the instantaneous vorticity obtained on the domain L/d = 10 at Reynolds numbers 
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Figure 3: Cylinder flow: time histories of the lift force on the cylinder at Reynolds numbers (a) Re = 60 
and (b) Re = 500. Results are obtained with Dq = in the open boundary condition. 


Reynolds number 

Domain 

Cd 

C'd 

Cl 

20 

T/d = 5 

2.294 

0 

0 


L/d= 10 

2.317 

0 

0 


L/d= 15 

2.317 

0 

0 


L/d = 20 

2.317 

0 

0 

100 

L/d = 5 

1.441 

8.491F;- 3 

0.261 


L/d= 10 

1.459 

7.631F;- 3 

0.254 


L/d= 15 

1.462 

7.700F;- 3 

0.253 


L/d = 20 

1.462 

7.714F;- 3 

0.253 


Table 1: Cylinder flow: effect of the domain size on the global flow parameters. Cd- drag coefficient or 
time-averaged mean drag coefficient; C^: rms drag coefficient; Cl- rms lift coefficient. 


Re = 30 and Re = 200, respectivefy. The results show a steady-state flow at Re = 30 and regular vortex 
sheddings at Re = 200. 

We have computed and monitored the forces acting on the circular cylinder. In Figure [3] we show a 
window of the time histories of the lift (i.e. the force component in the cross-flow y direction) at Reynolds 
numbers Re = 60 and 500. The force signals exhibit quite regular fluctuations about a zero mean value at 
these Reynolds numbers. 

Global flow parameters can be determined based on these force data. In Table [T] we have listed several 
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O Experiment (Roshko, Tritton, etc) 

O Dong & Karniadakis (2005), 3D, Re=10,000 
< Ma et al. (2000), 3D, Re=3,900 

□ Dong & Shen (2015), 2D 

A Current simulations, 2D 


O 



Re for 3D flow onset 


nf ^ .10’' ^ .10= ^ .lo" 

Re 


□ 

Dong & Shen (2015), 2D 


— Empirical relations, Norberg (2003) 

o 

Current simulations, 2D 


0.8 



(a) (b) 


Figure 4: Cylinder flow: Comparison of (a) drag coefficient and (b) rms lift coefficient versus Reynolds 
number between the current simulations, the experimental data, and the simulation results of 0. 


flow parameters for two Reynolds numbers Re = 20 and 100 obtained on different flow domains. They 
include: the drag coefficient Cd = t% 2 , where / is the time averaged drag (i.e. force component in x 
direction) and p = 1 is the fluid density; the root-mean-square (rms) drag coefficient Cj = where 

2 P^O 

f' 

is the rms drag; and the rms lift coefficient = i lr 2 , where /' is the rms lift. These data indicate that 
the size of the wake region has an influence on the simulation results, and that as the wake region becomes 
sufficiently large the flow parameters computed from the simulations remain essentially unchanged as the 
domain size further increases. It can be observed that the flow domain L/d = 10 is very close to the point 
where further increase in the domain size no longer results in significant changes in the results. In light 
of this observation, our subsequent discussions will be mainly based on the results obtained on the domain 
L/d =10. 


To demonstrate the accuracy of the method, we compare the force parameters computed from the current 
simulations with those from the experimental measurements and from other simulations in the literature. In 
FigurelHa) we plot the drag coefficient (Cd) as a function of Reynolds number from the current simulations. 


from a number of experiments 
simulations in 


11 , ^ 21 , 69 , 


61j . and from the simulations of 


m 


11 


■i 




Note that the 
are three- 


18| and in the current work are both two-dimensional, while those of 
dimensional. The current results correspond to the domain size L/d = 10 and Dq = in the open boundary 
condition 0. They agree with those of [l8| very well. Note that both the numerical algorithm and the 


outflow boundary condition in the current work are different from those of |18l |. In the two-dimensional 
regime the current results also agree well with the experimental data. But for Reynolds numbers where 
the physical flow has become three-dimensional (beyond about Re = 180), the current two-dimensional 
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Reference 

Re = 100 

Re = 200 

Braza et al. (1986) 

0.21 

0.55 

Engelman & Jamnia (1990) [1^ 

0.26 

- 

Meneghini & Bearman (1993) 

- 

0.54 

Beaudan & Moin (1994) [2] 

0.24 

- 

Zhang et al. (1995) [7^ 

0.25 

0.53 

Tang & Audry (1997) [66] 

0.21 

0.45 

Persillon & Braza (1998) 

0.27 

0.56 

Zhang & Dalton (1998) [7^ 

- 

0.48 

Kravchenko et al. (1999) 

0.22 

- 

Hwang & Lin (1992) [^ 

0.27 

0.42 

Franke et al. (1990) [2^ 

- 

0.46 

Karniadakis (1988) [3^ 

- 

0.48 

Newman & Karniadakis (1995) [^ 

- 

0.51 

Newman & Karniadakis (1996) [^ 

0.24 

- 

Dong & Shen (2010) [1^ 

~ 

0.501 

Dong & Shen (2015) [1^ 

0.254 

0.527 

Current simulation (domain L/d = 10) 

0.254 

0.526 

Current simulation (domain L/d = 20) 

0.253 

- 


Table 2: Cylinder flow: Comparison of rms lift coefficients at Re = 100 and Re = 200 between current 
simulations and other simulations from literature. 


simulations result in overly large drag coefficients compared to the experiments, and the discrepancy grows 
with increasing Reynolds number. 

Figure m^b) is a comparison of the rms lift coefficient Cl as a function of the Reynolds number between 
current simulations, the experiment of [^ . and the simulations of 


relation given by 


18|. The curves show the empirical 


53| based on several experimental sources, which exhibits a hysteresis around the Reynolds 


numbers where the two-dimensional to three-dimensional flow transition occurs. The lift coefficients from 


the current simulations and from 18| agree with each other almost exactly. In the two-dimensional regime 


the current results agree with the empirical relation from 


53j| reasonably well. In the three-dimensional 


regime, however, the current two-dimensional simulations grossly over-predict the lift coefficient, which is a 
well-known phenomenon about two-dimensional simulations (see e.g. |13l. Il5|b 

In Table [2] we have summarized the rms lift coefficients {Cl) for Reynolds numbers Re = 100 and 200 
from a number of two-dimensional simulations in the literature. We have also listed the Cl values on flow 
domains with L/d = 10 and 20 from the current simulations for comparison; see Table [1] for Cl values on 
the other domains at Re = 100. One can observe a spread in the Cl values from different simulations. The 
lift coefficients from the current work are well within the range of values from the literature. 

Let us next look into the effectiveness of the open boundary condition and the algorithm from Section 
[5] for dealing with the backflow instability. We will consider the cylinder flow at higher Reynolds numbers, 
ranging from Re = 2000 to Re = 10000. At these Reynolds numbers energetic vortices are observed to pass 
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Figure 5: Cylinder flow: snapshots of instantaneous velocity fields at Reynolds numbers (a) Re = 2000, 
(b) Re = 5000, and (c) Re = 10000. Velocity vectors are plotted on every fifth quadrature points in each 
direction within each element. 
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Figure 6: Time histories of the lift force on the cylinder at Reynolds numbers (a) Re = 2000, and (b) 
Re = 10000. 


through the outflow boundary and induce strong backflows in that region. This creates a severe instability 
issue, and makes the simulation immensely challenging. 

Thanks to the energy stability, the current open boundary condition provides an effective way for over¬ 
coming this instability. In Figure [5] we show distributions of the instantaneous velocity at three Reynolds 
numbers Re = 2000, 5000 and 10000. The results are obtained on the domain L/d = 10, with Dq = in 
the open boundary condition (|4]). Energetic vortices can be clearly observed at the open boundary (see e.g. 
Figure EKc)). 

We have performed long-time simulations at these Reynolds numbers using the current method. The 
long-term stability of the simulations is demonstrated by Figure |6l in which we show a window of the time 
histories (over 30 flow-through time) of the lift force on the cylinder at Reynolds numbers Re = 2000 and 
Re = 10000 obtained on the domain L/d = 10. At Re = 2000 the lift signal exhibits a modulation in its 
amplitude. At Re = 10000 the fluctuation becomes quite chaotic and the vortex-shedding frequency appears 
to vary over time. These results show that simulations using the current method are long-term stable in the 
presence of strong vortices and backflows at the outflow/open boundaries. 

In contrast, the boundary condition (|T0)) is observed to be unstable at these Reynolds numbers considered 


23 














(a) 



Figure 7: Time histories of the drag on the cylinder at Re = 10000 obtained using different Dq values in 
OBC: (a) Dq = 0, and (b) Dq = 

here. The computation instantly blows up as the vortices hit the outflow/open boundary. 

Let us next look into the effects of the Do value in the open boundary condition (jH) on the simulation 
results. By an analogy between the current open boundary condition and the usual convective boundary 
condition (HU, we observe that should correspond to a convection velocity scale Uc at the outflow 
boundary, i.e. = Uc, as is discussed in Section [2j The simulation results for the cylinder flow presented 
so far are obtained with a value -^ = Uc = Uq, where Uq is the average velocity at the outflow boundary. 

We have observed that the variation in the Dq value in the open boundary condition (|4]) has little or 
essentially no effect on the global physical quantities such as the forces on the cylinder. This is illustrated 
in Figure [7] with the time histories of the drag at Reynolds number Re = 10000 corresponding to two 
values DqUo = 0 (or Uc = oo) and DqUo = 1 (or Uc = Uq) in the open boundary condition (jd]). The two 
drag signals exhibit qualitatively similar characteristics. A quantitative comparison of the global physical 
quantities corresponding to several Dq values is given in Table [3l The table includes the data for the mean 
drag coefficient {Cd), rms drag coefficient (C^), and the rms lift coefficient (Cl) for Reynolds numbers 
Re = 20, 100 and 10000 on the flow domain with L = lOd. We have considered several Dq values in these 
tests, corresponding to DqUq = 0, 0.5, 1, 2 and 5, or equivalently Uc = oo, 2Uo, Uq, Hf- and One can 
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Reynolds number 

DqUq 

Cd 

C'd 

Cl 

20 

0.0 

2.317 

0.0 

0.0 


0.5 

2.317 

0.0 

0.0 


1.0 

2.317 

0.0 

0.0 


2.0 

2.317 

0.0 

0.0 


5.0 

2.317 

0.0 

0.0 

100 

0.0 

1.459 

7-627E - 3 

0.254 


0.5 

1.459 

7.630D - 3 

0.254 


1.0 

1.459 

7.631D- 3 

0.254 


2.0 

1.459 

7.639D - 3 

0.254 


5.0 

1.459 

7.656D - 3 

0.254 

10000 

0.0 

1.862 

0.449 

1.483 


0.5 

1.881 

0.421 

1.506 


1.0 

1.843 

0.442 

1.474 


2.0 

1.893 

0.432 

1.504 


5.0 

1.858 

0.446 

1.474 


Table 3: Effect of Do in OBC on global ffow parameters: drag and lift coefficients for cylinder flow obtained 
with different Dq values. Cd- drag coefficient or time-averaged mean drag coefficient; C'^: root-mean-square 
(rms) drag coefficient; Cl- rms lift coefficient. 


observe that at Re = 20 and Re = 100 these global physical quantities obtained using these several Dq 
values are exactly or almost exactly the same. At Re = 10000, they are also close for different Dq values, 
with a maximum difference of about 2.7% for Cd, about 6.7% for C^, and about 2.1% for Cl- These results 
suggest that the value of Dq in the open boundary condition (|T|) has little effect quantitatively on the physical 
quantities of the flow. 

The main effect of Dq appears to be on the qualitative features of the flow, such as the smoothness of the 
velocity field, in regions local to the outflow boundary. We observe that the current open boundary condition 
([^ with appropriate Dq > 0 result in smoother velocity distributions at the outflow boundary and can allow 
the vortices to cross the outflow boundary more smoothly and in a more natural way when compared to the 
boundary condition with Dq = 0, which corresponds exactly to the boundary condition OBC-C studied in 


18| . These observations suggest that one should choose in accordance with the convection velocity scale 


of the vortices at the outflow boundary, which for the current problem is close to Uq, in order to improve 
the qualitative characteristics of the computed flow near the outflow boundary. A larger deviation of the 
chosen value in the open boundary condition from the actual vortex-convection velocity at the outflow 
boundary may lead to qualitatively poorer flow patterns near the outflow boundary. 

3.3 Jet in Open Domain 


In this subsection we apply the current method to simulate a jet in an open domain. This is a type of flow 
different than that of the previous subsection. The open domain boundaries combined with the physical 
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Figure 8: Configuration of jet in an open domain. 


instability of the jet at large Reynolds numbers make this type of flows very challenging to simulate. We 
will again consider two-dimensional simulations. 

Figure [5] is a sketch of the flow configuration. We consider a jet with an inlet diameter d, which is 
contained in a rectangular domain, —2.5(i ^ a; ^ 2.bd and 0 ^ 1 / T.Sd. The bottom of the domain (y = 0) 
is a solid wall, and the jet inlet is located in the middle of the wall, covering —0.5d ^ x ^ 0.5d. The other 
three boundaries (top, left and right) of the domain are open, where the fluid can freely leave or enter the 
domain. We assume that the jet velocity at the inlet has the following profile 



{H{x, 0) — H{x, d/2)) tanh —+ {H{x, —d/2) — H{x, 0)) tanh 


(41) 


V2e ^ ^ ^ ^ V2e \ 

where {u,v) are components of the velocity u in x and y directions, Uq is a velocity scale, and e = ^. 
H{x,a) is the Heaviside step function, assuming a unit value if x ^ a and vanishing otherwise. Note that 
with this profile the inlet flow has a bulk velocity Uq. We assume that there is no external body force for 
this problem. 

We normalize all the length variables by the jet inlet diameter d, and all the velocity variables by Uq. 
The Reynolds number is therefore given by equation (1401) . in which d and Uq have meanings particular to 
this problem. 

The domain has been discretized using 600 equal-sized quadrilateral spectral elements, with 20 elements 
in the x direction and 30 elements in the y direction. We impose the velocity Dirichlet condition ([3]) on the 
bottom side of the domain, with the boundary velocity w(x, f) = 0 in the wall region and set according to 
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Figure 9: Jet in open domain: snapshots of instantaneous velocity fields at (a) Re = 2000, (b) Re = 5000, 
and (c) Re = 10000. Velocity vectors are plotted on every nineth quadrature points in each direction within 
each element. Results are obtained with Dq = in the open boundary condition. 


equation (1411) in the jet inlet. On the top, left and right sides of the domain we impose the open boundary 
condition (|4|) with ff, = 0 and S = The majority of simulations in this subsection are performed using 
Dq = in the open boundary condition. Several other Dq values have also been considered in selected 
cases for comparison. 

We integrate the Navier-Stokes equations (lla|) - (llbl) in time using the algorithm described in Section [2j 
Long-time simulations have been performed at three Reynolds numbers: Re = 2000, 5000 and 10000. In the 
simulations we employ an element order 12 for each element at the two lower Reynolds numbers, and an 
element order 16 for each element at Re = 10000. The non-dimensional time step size is = 2.5 x 10“'* 
for Re = 2000 and 5000, and = 2.0 x 10“* for Re = 10000. 

The flow characteristics are illustrated by snapshots of the instantaneous velocity shown in Figure IHl at 
these Reynolds numbers. These results are obtained using Dq = in the open boundary condition (jT]) . At 
Re = 2000, the jet appears to be stable within a distance of at least 3d downstream of the inlet (FigurelHl^a)). 
Beyond this region, the jet becomes unstable and a pair of vortices forms, which travels downstream along 
with the jet and eventually crosses the upper open boundary and exits the domain. The velocity field appears 
symmetric with respect to the jet centerline at this Reynolds number. As the Reynolds number increases to 
Re = 5000, the stable region immediately downstream of the inlet becomes shorter (approximately 2d, see 
FigureinKb)), and the velocity distribution has lost the symmetry with respect to the jet centerline. Pairs of 
vortices can be observed to form, wrapping around the jet at different downstream locations. At Re = 10000, 
The stable region downstream of the inlet becomes evern shorter (about d), and the vortices forming along 
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Figure 10: Jet in open domain: Time histories of the vertical force component exerting on the wall at 
Reynolds numbers (a) Re = 2000, (b) Re = 5000, and (c) Re = 10000. Results correspond to Dq = -^ in 
the open boundary condition. 


the jet are notably more numerous (Figure [HIc)). 

The jet simulations using the current method are long-term stable, even in the presence of backflows 
or vortices (see e.g. Figure lljc)) at the open domain boundaries. This is demonstrated by Figure HOl in 
which we show a window of the time histories of the vertical force acting on the bottom wall for these three 
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Figure 11: Jet in open domain {Re = 5000): time histories of the vertical force component obtained with 
different Uq values in OBC: (a) DqUo = 0 and (b) DqUo = 2.0. These can be compared with Figure [TUT bl. 
which corresponds to DqUo = 1.0 at the same Reynolds number. 


Reynolds numbers. Note that the horizontal force {x component) on the wall is essentially zero. These results 
again correspond to Dq = in the open boundary condition ((T|) . One can observe that the force signals are 
highly unsteady, and the fluctuations appear more energetic and chaotic with increasing Reynolds numbers 
fFigure irUT cll. The large temporal window shown here, over 500-^ or approximately 67 flow through times, 
signihes the long-term stability of our simulations. 

The term containing 0o(n, u) in the open boundary condition (|4|) is critical to the stability of the current 
simulations. For comparison, we have also employed the boundary condition (HOD on the open boundaries 
for simulations of the current problem, and observe that the computation is unstable for all the Reynolds 
numbers considered here. The simulations blow up instantly as the vortices reach the open boundaries. 

The results shown so far have been obtained using Dq = -^ in the open boundary condition (|4|) . We next 
concentrate on the Reynolds number Re = 5000, and look into the effect of Dq on the simulation results. 
We have considered several Dq values in the open boundary condition: DqUq = 0, 0.5, 1.0, and 2.0. Figure 
m shows time histories of the vertical force on the wall at Re = 5000 obtained using Dq = 0 and Dq = 
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DoUo 


fy 


f' 

•fy 


0.0 

-1.621F;- 6 

2MIE- 

3 

QMbE- 

3 

0.5 

1.203F;- 6 

2MAE - 

3 

5.713F;- 

3 

1.0 

-9.599F;- 7 

2.729E- 

3 

5.549F;- 

3 

2.0 

2.829E - 7 

2.910F;- 

3 

5.541F;- 

3 


Table 4: Mean and rms forces on the wall for the jet problem at Re = 5000 obtained using different 
Dq values in OBC. f ^ and fy denote the time-averaged mean forces in horizontal and vertical directions, 
respectively, fy is the rms force in the vertical direction. 


in the open boundary condition (|4]). These plots can be compared with Figure [TOT b'). which corresponds to 
Dq = -jh- at the same Reynolds number. These force signals appear qualitatively similar. 

To provide a quantitative comparison, we have listed in Table 0] the time-averaged mean and rms forces 
acting on the wall corresponding to different Dq values in the open boundary condition. The mean force in 
the horizontal direction (f^) is essentially zero for all cases. Both the mean and rms forces in the vertical 
direction obtained using different Dq values are close, with a maximum difference of about 6.6% for the mean 
vertical force and a maximum difference of about 9.5% for the rms vertical force. The Dq effect on the global 
flow quantities (such as forces) observed here for the jet flow is approximately in line with the observations 
for the cylinder flow in Section [321 But the differences observed for this problem appear slightly larger. 


4 Concluding Remarks 


The main contributions of the current paper are two-fold: 


• We have presented a new type of energy-stable open boundary condition for incompressible flow sim¬ 
ulations. This boundary condition ensures the energy stability of the system, and in some sense can 
be analogized to the usual convective boundary condition. The current open boundary condition can 
be reduced to that of [l^ if the inertia term involved herein vanishes (by setting Dq = 0). Note that 
the current boundary condition can be generalized to a famil y o f convective-like energy-stable open 
boundary conditions with an idea similar to that employed in 


18|; see the remarks in Section [2T 


• We have presented an efficient numerical algorithm for the proposed open boundary condition. The 
key issue here lies in the numerical treatment of the inertia term involved in the current open bound¬ 
ary condition. Our algorithm combines a rotational velocity-correction type splitting strategy for the 
Navier-Stokes equations with two different implicit approximations of the inertia term in the open 
boundary condition for the pressure and velocity sub-steps. The algorithm leads to Robin-type condi¬ 
tions for both the discrete pressure and the discrete velocity on the outflow/open boundary. 
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We would like to point out that the open boundary condition proposed herein can be considered as the 
energy-stable version of a combined traction-free and convective boundary condition (see equation OT). 
As discussed in Section [2A] this condition will be reduced to the traction-free condition if the inertia term 
vanishes (i.e. Dq = 0)i £md if p = 0 is imposed on the outflow boundary it will be reduced essentially to the 
usual convective boundary condition. 

The method developed in the current work for dealing with outflow/open boundaries exhibits favorable 
properties when compared with those of [l4, 1^ in at least two aspects. First, it can lead to qualitatively 
smoother flow patterns at or near the outflow boundary, because the current boundary condition (when 
Dq > 0) provides a control over the velocity on the outflow boundary. This has been demonstrated by the 
numerical simulations of Section [3] Second, it provides a simpler implementation (with Dq > 0) with 
spectral-element and finite-element type methods. This is because the current method essentially imposes a 
discrete Robin-type condition for both the velocity and the pressure on the outflow boundary, and requires 
only an update to the coefficient matrix by a surface integral in the implementation. In contrast, with the 
methods of [l4, 18| a pressure Dirichlet type condition is imposed on the outflow boundary, and a projection 
to the H^{dilo) space of the pressure Dirichlet data will be required with elements in the implementation. 
This projection is more involved than the evaluation of the surface integral required by the current method. 

The effectiveness of the current method has been demonstrated by extensive numerical simulations. 
Comparison with the experimental data shows the accuracy of the current method. At higher Reynolds 
numbers when backflows and strong vortices occur at the outflow/open boundaries, numerical results have 
demonstrated the long-term stability using the current method. It is observed that the method allows the 
vortices to discharge from the domain in a fairly natural fashion, even at quite high Reynolds numbers (up 
to Re = 10000 tested here). 

We anticipate that the method developed herein will be instrumental in numerical studies of wakes, jets, 
shear layers, and other types of flows involving physically unbounded domains, especially for high Reynolds 
numbers. It would be very interesting to extend the idea and develop analogous boundary conditions for 
moving domains such as in an arbitrary Lagrangian Eulerian context. This would be important and useful to 
applications such as vortex/flow induced vibrations. Future research will address such and related problems. 
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